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Analysis  of  1-D  Moment  Equations  for  Immiscible  Flow 

Kenneth  D.  Jarman  and  Thomas  F.  Russell 

Abstract.  We  derive  and  analytically  and  numerically  solve  statistical  mo¬ 
ment  equations  for  immiscible  flow  in  porous  media  in  the  limit  of  zero  capil¬ 
lary  pressure,  with  application  to  secondary  oil  recovery.  Mean  and  variance  of 
(water)  saturation  exhibit  a  bimodal  character;  two  shocks  replace  the  single 
shock  front  evident  in  the  classical  Buckley-Leverett  saturation  profile. 


1.  Introduction 

Subsurface  geologic  properties  at  field  scales  are  uncertain,  and  are  often  de¬ 
scribed  statistically  in  practice.  Flow  profiles  in  such  porous  media  are  uncertain, 
and  statistical  flow  outcomes  are  appropriate.  We  are  primarily  interested  in  mean 
behavior  and  a  measure  of  the  uncertainty  about  this  mean. 

A  “zeroth-order”  model  of  mean  flow  with  averages  of  geologic  properties  ig¬ 
nores  correlations  between  flow  variables.  Monte  Carlo  simulations  of  many  realiza¬ 
tions  of  geologic  properties  to  estimate  moments  requires  much  computation  time 
and  careful  sampling  techniques  [8],  [9],  [28].  Macrodispersion  theories  in  contam¬ 
inant  transport  capture  a  first-order  effect  of  fluctuations  via  covariance  functions 
in  PDFs  for  the  mean  concentration  [5],  [10]. 

We  derive  second-order  PDFs  for  the  covariance  functions  and  the  mean  flow, 
and  solve  for  these  moments  simultaneously.  The  fundamental  problem  of  closure 
of  the  system  is  addressed  by  a  perturbation  argument.  The  resulting  moment 
equations  directly  approximate  the  local  mean  and  covariance  functions,  for  general 
boundary  conditions  and  general  stochastic  geology  [31]. 

1.1.  Applications.  A  statistical  description  of  subsurface  flow  is  of  particular 
interest  for  secondary  oil  recovery.  The  principal  difficulty  is  a  non-convex  nonlin¬ 
ear  flux  function  in  an  advection  equation  that  leads  to  discontinuous  solutions. 
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Standard  pressure-saturation  equations  for  2-D  horizontal  flow  of  two  immiscible 
fluids  in  porous  media,  in  the  limit  of  vanishing  capillary  pressure  {jpc  =  0) ,  are 

(1.1)  ^  v(x)  = -K(x)  V/i(x),  V  •  v(x)  =  0, 

(1.2)  5ts(x,  t)  +  V  •  [/(s(x,  t))v(x)]  =  0. 

These  are  taken  to  be  valid  from  laboratory  (centimeters)  to  field  scales  of  reservoir 
depth  (10-100  meters)  and  length  (100-1000  meters).  Hydraulic  conductivity  K 
may  be  an  anisotropic  tensor;  here,  for  simplicity,  it  will  be  an  isotropic  scalar  K. 
Assume  also  that  K  depends  weakly  on  (water)  saturation  s  [1],  [3].  Apply  (1.1)- 

(1.2)  to  the  flow  of  oil  and  water,  for  arbitrary  fluid  mobilities.  Denote  the  total 
velocity,  a  scaled  total  volumetric  flux  of  both  fluids,  by  v,  hydraulic  water  head 
by  h,  and  porosity  (assumed  constant)  by  (p.  The  fractional  flow  function  f{s), 
for  Pc  =  0,  represents  the  fraction  of  v  due  to  water.  It  is  typically  S-shaped;  we 
use  here  a  functional  form  arising  from  quadratic  relative  permeabilities  (see  [1]), 
though  our  method  does  not  depend  on  any  such  specific  choice. 

Capillary  pressure  regularizes  sharp  fronts  caused  by  the  nonlinear  advection 
term.  To  obtain  a  linear  approximation  to  this  effect,  add  e£)V^s(x,  t)  with  eu  >  0 
to  the  right  side  of  (1.2).  Letting  eu  — >■  0  defines  the  vanishing-viscosity  solution 
[25],  which  is  the  one  we  seek. 

As  is  standard  in  subsurface  applications,  let  Y  =  Ini^  be  a  random  held 
with  prescribed  mean  and  covariance  functions;  e.g.,  it  is  often  claimed  that  Y  is 
multivariate  Gaussian,  based  on  empirical  observations  [6]  (our  method  does  not 
depend  on  this).  Through  (1.1)-(1.2),  v  and  s  are  thus  random  fields.  No  other 
underlying  sources  of  uncertainty  are  considered  in  this  study. 

Under  the  assumptions  stated  above,  with  steady  boundary  conditions,  a  steady 
V  can  be  determined  from  (1.1).  We  evolve  s  from  the  stochastic  PDE  (1.2), 
assuming  v  is  known.  Moments  of  v  and  h  can  be  estimated  from  established 
theory  ([29]  and  [31]  use  moment  equations).  We  seek  to  combine  analytical  and 
numerical  techniques  to  model  the  propagation  of  uncertainty  from  an  underlying 
random  held  y(x),  through  v(x),  to  the  solution  s{x,t). 

1.2.  Previous  work.  Existing  work  on  moment  differential  equations  (MDEs) 
focuses  mostly  on  advection  equations  with  linear  flux  functions  [1 1] ,  and  some  non¬ 
linear  subsurface  flow  equations  of  a  form  different  from  (1.2)  [2],  [27],  [29],  [30]. 

Langlo  and  Espedal  [16],  [17]  presented  a  macrodispersion  approach  for  the 
stochastic  version  of  (1.2).  The  flux  function  is  expanded  in  a  Taylor  series,  and 
high-order  terms  are  neglected;  then  standard  techniques  represent  macrodispersiv- 
ity  as  a  function  of  flow  velocity  covariance.  Zhang,  Tchelepi,  and  Li  take  advantage 
of  the  steady  velocity  held,  and  transform  2-D  flow  to  1-D  Lagrangian  flow  along 
streamlines  [32],  [33].  Then  they  formulate  integral  equations  for  moments  from 
ensemble  averages  over  the  streamtubes. 

An  Eulerian  MDE  approach  has  been  successful  for  single-  and  multiphase 
pressure  and  velocity  equations  [29],  and  a  natural  next  step  is  to  extend  the 
theory  from  flux  equations  to  transport  equations.  This  framework  differs  from 
streamtubes  not  only  in  formulation,  but  also  in  that  the  MDEs  need  no  velocity- 
distribution  assumption,  and  an  extension  to  transient  velocity  fields  is  relatively 
straightforward.  The  approach  applies  to  any  probability  distribution  of  geologic 
properties  and  any  correlation  function,  and  does  not  require  stationarity.  Other 
stochastic  theories  generally  require  such  restrictions. 
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Equations  are  derived  in  §2.  In  §3,  we  reduce  the  mean  and  variance  equations 
to  a  simple  form.  These  are  shown  to  be  strictly  hyperbolic  for  1-D  flow,  and  an 
analytical  solution  is  given  in  §4.  Classification  is  briefly  discussed  for  2-D  equa¬ 
tions.  Details  of  the  results  outlined  here  are  presented  in  [13].  We  conclude  with 
an  evaluation  of  our  results  and  their  practical  implications,  and  a  brief  overview 
of  additional  questions  that  warrant  future  investigation. 

2.  Moment  Equations 

We  compare  two  different  approaches  for  statistical  MDEs  of  1-D  flow.  In  §2.1, 
we  expand  fields  in  the  form  u  =  {u)  +  5u  and  the  resulting  MDEs  are  closed  by 
neglecting  products  of  ‘d’  terms.  Analogous  equations  in  §2.3  result  from  a  full 
asymptotic  expansion.  We  show  in  §4  that  the  latter  yields  moments  that  violate 
physical  constraints.  In  both  cases,  random  fluctuations  in  E  =  Ini^  are  assumed 
small:  ay  <C  1.  We  can  immediately  generalize  to  higher  dimensions  via  vector 
notation.  Moments  of  h  and  v  are  assumed  known  from  (1.1)  and  moments  of  Y. 
The  1-D  saturation  equation  (1.2),  with  initial  data  s(x,0)  =  g{x),  is 

(2.1)  dts  +  d^{f{s)v)  =0, 

We  assume  that  g  is  known  with  certainty.  Solutions  are  defined  in  terms  of  van¬ 
ishing  viscosity  as  in  §1.1;  henceforth,  this  is  tacitly  understood. 

2.1.  Two-term  expansion.  Moment  equations  may  be  derived  in  a  number 
of  ways.  Eor  examples  of  commonly  used  methods  applied  to  various  models  in 
subsurface  flow  and  transport,  see  [5],  [7],  [10],  [11],  [24],  [29].  Here  we  apply  a 
standard  approach,  separating  mean  fields  from  random  fluctuations. 

Let  (•)  denote  the  expectation  operator,  defined  by 

(2.2)  (V>)  =  /  V’(w)  dP(w) 

Jn 

for  any  integrable  function  ip  :  ^  R  on  the  sample  space  D  with  probability 

measure  P.  We  omit  reference  to  w  in  what  follows. 

The  random  field  Y  is  decomposed  into  deterministic  mean  plus  random  fluc¬ 
tuation:  Y  =  (Y)  +  6Y.  Each  field  dependent  on  Y  is  represented  similarly: 

(2.3)  h{x)  =  {h)  (x)  +  5h{x),  v{x)  =  {v)  (x)  -|-  5v{x), 

s{x,t)  =  (s)  {x,t)  +  6s{x,t). 

Recall  that  we  only  need  the  decompositions  of  v  and  s  here.  Next,  the  fractional 
flow  function  is  expanded  in  a  Taylor  series  around  {s): 

(2.4)  /(.)  =  /((.))  +  n(s))Ss  +  \n(s))Ss^  +  •  •  • 

So  far,  we  make  no  assumption  regarding  the  size  of  5s  relative  to  (s). 

To  obtain  the  mean-saturation  equation,  apply  the  operator  (•)  to  (2.1): 

(2.5)  dt  {s)  +  d,  [/((.))  (u)  +  /'((.))  (d.  du)  +  1/" ((.))  (u)  (d.^) 

+  (Ss^Sv)  +  ^rUs))  {v)  {6s^)  +  •  •  •  ]  =  0. 
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The  flux  terms  include  a  nonlinear  advective  mean,  two  covariances,  and  higher- 
order  moments.  For  the  saturation-fluctuation  equation,  subtract  (2.5)  from  (2.1): 

dt5s  -h  5a,  |^/((s))du  -h  f{{s))5s  (v)  +  f'{{s))  {6s  6v  -  {6s  6v)) 

(2.6)  +^/"((s))  {v)  {6s^  -  (6s^})  +  ^f"((s))  (6s^6v  -  (6s^6v})  -h  •  •  •]  =  0. 

We  derive  equations  for  the  unknown  covariance  functions  {6s6v)  and  {6s^),  using 
the  following  additional  notation.  The  independent  variables  are  x  and  t  except 
where  noted,  and  ■\y  denotes  the  replacement  of  x  by  some  y  different  from  x. 
It  is  convenient,  and  useful,  to  derive  equations  for  the  more  general  two-point 
covariances  and  rather  than  for  one-point  covariances. 

To  obtain  an  equation  for  the  saturation-velocity  covariance  (dsdu|j^),  multiply 

(2.6)  by  6v{y)  and  apply  (•).  This  results  in: 

dt  +  5a,  [/((s))  {6v  6v\y)  +  f'{(s))  (v)  (6s  6v\y)  +  f'{(s))  (6s6v  6v\y) 

(2.7)  +^/"((s))  {v)  +  ^/"((s))  {6s‘^6v  -h  •  •  •  ]  =0. 

Similarly,  multiplying  (2.6)  by  6s{y,t),  and  using  the  identity^  6sdt6s\y  -\- 
6s\ydt6s  =  5t(dsds|j^),  yields  this  equation  for  the  two-point  saturation  covariance: 

dt  (dsdsly)  +  5a,  [/((s))  +  f{{s))  {v)  {6s  dslj,)  -h  f'{{s))  {6s  6v 

(2.8)  +\n{s))  {v)  {6s^6s\y)  +  i/" ((.))  d.|,)  +  •  •  •  ] 

+53/[/((«l3/»  dulj,)  -|-/'((s|3,))  (ulj,)  (dslj^ds)  -|-/'((s|3,))  {{6s  du)|3,ds) 

+  ln{s\y))  {v\y)  (ds^l.ds)  +  \f'{{s\y))  {{6s^6v)\y6s)  +  •••]=  0. 

2.2.  Closure  by  perturbation  argument.  If  ay  <C  1  so  that  fluctuations 
and  their  derivatives  may  be  assumed  small  relative  to  the  means,  and  if  /  is 
smooth,  then  we  can  approximate  (2.5)-(2.8)  by  a  closed,  coupled  system.  Defin¬ 
ing  Cs^{x,y,t)  =  (dsdulj,),  =  (dsdslj,) ,  =  (dudulj,) ,  (s)  Ij,  =  {s){y,t),  and 

Csy{x,  y,  t)  =  Csy{y,  X,  t),  the  resulting  system  is 

(2.9a)  dt  (s)  +  5a,  [f{(s))  (v)  +  f'{(s))a,^  +  ^f"{(s))a‘^,  (v)  ]  =  0, 

(2.9b)  dtCsv  +  5a,  [/((s))Ca,  -h  f'{{s))  (v)  Csa,]  =  0, 

(2.9c)  dtCs  +  5a,  ^f{{s))Csv  +  /'((«))  {v)  Cs]  -h  dy  ^f{{s'^y)Csv  +  f{{s)\y)  (uIj,)  Cs]  =  0. 

Initial  data  are  (s)  (x,  0)  =  g{x),  Csv{x,y,0)  =  Cs{x,y,0)  =  0;  recall  that  {v)  and 
{6v6v\y)  are  assumed  known.  Both  (2.9b)  and  (2.9c)  have  advective  flux  terms,  are 
coupled  to  the  mean  equation  (2.9a),  and  are  first-order  in  ay.  This  is  consistent 
with  the  approximation  to  (2.9a),  which  is  second-order  in  ay. 

Another  common  closure  argument,  which  may  be  called  a  Gaussian  assump¬ 
tion,  might  be  applied  here.  For  example,  for  the  linear  case  {f{s)  =  s),  (2.9)  is 
exact  if  one  assumes  that  velocity  and  saturation  are  jointly  multivariate  normal  [7]. 

^The  identity  is  not  valid  in  a  strong  sense  for  discontinuous  solutions.  Recall,  however,  that 
we  define  solutions  in  terms  of  the  (smooth)  viscous  solution,  in  the  limit  en  0. 
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This  follows  without  a  perturbation  argument.  The  linear  case  under  this  assump¬ 
tion  is  studied  in  [11],  where  the  effect  of  small-scale  diffusion  is  included.  Using 
simulations,  we  found  that  the  Gaussian  assumption,  even  for  transformations  of 
saturation  and  velocity  fields,  is  inappropriate  for  this  problem  [12]. 

Note  that  the  mean  equation  (2.9a)  contains  the  functions  Gsv  =  Csv{x,  x,  t)  and 
CTj  =  Cs{x,x,t)  rather  than  Cs^{x,y,t)  and  Cs{x,y,t).  This  mix  of  one-point  and 
two-point  covariance  functions  prevents  us  from  immediately  treating  the  MDEs  as 
classically  hyperbolic,  even  though  they  can  be  put  in  conservation-law  form.  Also, 
independent  variables  x  and  y  are  permuted  in  (s)  and  in  (2.9c).  In  general, 
(s)  \y  ^  (s),  and  Cgy  ^  Csv  We  address  these  issues  in  §3. 

We  refer  to  this  problem  as  “l-D,”  even  though  the  covariance  functions  involve 
two  spatial  variables,  and  the  saturation  covariance  has  fluxes  in  both  directions. 
The  two  variables  represent  two  different  points  in  the  same  1-D  domain.  Similarly, 
the  “2-D”  problem  has  four  spatial  coordinates. 

2.3.  Infinite  expansion.  In  this  alternative  derivation,  flow  variables  head 
/i(x),  velocity  v{x)  and  saturation  s{x,t)  are  represented  by  formal  infinite  pertur¬ 
bation  series  expansions  in  powers  of  a  parameter  e: 

OO  OO  OO 

(2.10)  /i  =  ^e"/i„(x),  i;  =  ^e"u„(x),  s  =  ^  e"s„(x,  t). 

n=0  n=0  n=0 

The  expansion  parameter  e  =  ay  is  shown  to  be  appropriate  within  the  context  of 
flow  equations  (1.1)  [5,  pp.  184-190],  [31].  For  example,  for  single  phase,  stationary 
uniform  mean  flow  in  1-D,  a^  is  approximated  by  =  v^ay- 

Moment  equations  analogous  to  (2.9)  can  be  derived  [13]: 

(2.11a)  dtso  +  [/(so)t^o]  =  0,  dt  (si)  -h  [uo/'(so)  (si)  ]  =  0, 

(2.11b)  dt  (S2)  +  dx  |^/(so)  (^^2)  +  f'{so)(^sv  +  f'{so)  (S2)  vq  +  =  0, 

(2.11c)  dtCsv  +  dx  [/(so)g  +  /'(so)^^oCs^;]  =  0, 

(2. lid)  dtCs  +  dx  [fiso)csv  +  /'(so)^^oCs]  +  dy  [/(  ^{)\y)^sv  +  /'(sol3/)^^ol3/Csj  —  0. 

Initial  data  are  given  by  So(a^!0)  =  g{x),  (si)(x,0)  =  (s2)(a^,0)  =  Csy{x,y,Q)  = 
Cs{x,  y,  0)  =  0.  The  second-order  mean  is  sq  +  e  (^i)  +  {s2)-  Note  that  the  argu¬ 

ment  of  is  so!  the  zeroth-order  mean.  The  system  is  in  fact  closed  again  using 
a  perturbation  argument,  but  now  this  argument  is  contained  in  the  assumption 
that  the  formal  power  series  in  e  converges.  Thus,  throughout  §2,  second-order 
equations  are  closed  by  assuming  that  heterogeneity  is  weak  (ay  <C  1). 

3.  Classification  of  reduced  equations 

We  exploit  special  structure  in  the  1-D  equations  that  allows  simplification  and 
classification  of  (2.9)  and  (2.11).  For  classification,  we  need  the  following 

Definition  3.1.  Let  u{x,y,t)  :  x  ]i_,_  ->•  M",  F  :  M"  ->•  M",  and  G  :  M"  ->• 

ffi".  A  system  of  equations  in  conservation-law  form  is  given  by 

(3.1)  dju -h  5a,F(u) -h  53,G(u)  =  0,  u(x,  0)  =  uo(x). 

This  system  is  hyperbolic  if  the  eigenvalues  of  ci  ilF(u)  -|-  C2  DG{m)  are  real  for  all 
Cl,  C2  G  ffi,  where  DF  and  DG  are  the  Jacobian  matrices  of  F  and  G  [20]. 
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We  expect  the  equations  to  be  nearly  hyperbolic  since  the  original  deterministic 
PDE  (1.2)  is  hyperbolic,  but  the  systems  (2.9)  or  (2.11)  as  given  cannot  even 
be  written  in  conservation- law  form  due  to  the  inconsistencies  mentioned  above. 
However,  we  show  that  both  sets  of  moment  equations  reduce  to  hyperbolic  systems 
of  conservation  laws  and  yield  analytic  solutions.  The  key  discovery  that  allows  this 
reduction  is  a  relationship  between  covariance  functions  in  1-D. 

Observe  that  velocity  v  is  constant  in  1-D,  so  it  is  merely  a  random  variable 
rather  than  a  random  field.  Thus,  velocities  at  any  two  points  in  space  are  perfectly 
correlated  (they  are  the  same  random  variable).  Equivalently,  the  velocity  corre¬ 
lation  length  is  infinite.  This  infinite  correlation  length  characterizes  the  principal 
difference  between  stochastic  subsurface  flow  in  one  and  two  space  dimensions  [12]. 

3.1.  Two-term  expansion.  The  MDEs  (2.9)  reduce  to  a  system  of  hyper¬ 
bolic  PDEs.  The  first  step  to  this  end  addresses  the  inconsistencies  mentioned 
earlier. 

Both  expansion  terms  {v)  and  dv  are  constant:  by  applying  the  expectation 
operator  to  d^v  =  0  we  obtain  {v)  =  0,  so  that 


0  =  dxV  =  dx  {v)  +  dxSv  =>  dxSv  =  0. 

This  implies  that  Cy{x,  y)  is  also  constant,  and  that  Csy{x,  y,  t)  is  independent  of  its 
second  argument.  Consequently,  is  identical  to  the  second-order  approximation 
to  velocity  variance  and  Csv{x,y,t)  is  identical  to  asv{x,t). 

This  last  identity  removes  the  inconsistency  of  having  asv  instead  of  Cgv  We 
still  have  in  (2.9a),  instead  of  c^,  and  we  have  {s)  \y  and  Cgy  in  (2.9c).  A  key 
variance-covariance  relationship,  follows:  in  1-D,  the  saturation  profile 

is  completely  determined  by  knowledge  of  the  velocity,  for  any  positive  time.  Thus, 
saturation  and  velocity  are  perfectly  correlated. 

We  divide  (2.9b)  by  >  0,  and  retain  only  the  first  two  equations  in  (2.9) 
in  the  following.  Replace  by  and  by  to  reduce  the  system  to  the 

following  new  equations  for  mean  and  standard  deviation  of  saturation: 


(3.2)  St 


f{{s))  +  crvf{{s))crs  +  ^  {v)  f'{{s))c 
^vf{{s))  +  {v)f'i{s))as 


Dependence  on  the  second  space  variable  y  has  been  eliminated.  Thus,  (3.2)  is  in 
conservation-law  form,  with  u{x,t)  =  {{s)  ,<Ts),  and  flux  function 


_  [  (v)  /((^))  +  +  i  (u)  /"((.)). 


Lemma  3.2.  The  moment  equations  (2.9)  reduee  to  a  hyperbolie  system. 
Proof.  It  remains  to  show  that  (3.2)  is  strictly  hyperbolic.  The  Jacobian 

matrix  DF((s) ,  CTs)  =  (  |  has  entries 

V  Ji2  J22  J 

111  =  {v)  f'{{s))+(Jvf''{{s))<Js+'^{v) 

112  =  f^vf'i(s))  +  (v)  f''i(s))cj,,  i22  =  (v)  f'i(s)). 

DF  is  symmetric;  thus  the  eigenvalues  of  DF  are  real  [26],  hence  (3.2)  is  hyper¬ 
bolic.  In  fact,  we  show  that  (3.2)  is  strictly  hyperbolic;  i.e.,  DF  has  a  complete 
set  of  linearly  independent  eigenvectors  [20].  If  the  eigenvalues  are  distinct,  then 
they  must  correspond  to  independent  eigenvectors.  Thus,  we  may  assume  that  DF 
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has  only  one  distinct  eigenvalue.  This  occurs  when  the  discriminant  of  the  char¬ 
acteristic  polynomial  vanishes:  {ju  —  ^22)^  +  4:ji2  =  0.  But  this  would  imply  that 
ji2  =  0,  hence  DF  is  a  multiple  of  the  identity,  and  it  again  has  a  complete  set  of 
eigenvectors.  The  conclusion  follows.  □ 

In  addition,  the  eigenvalues  are  distinct  unless  {{s)  ,as)  G  {(0,0),  (1,0)}  [12]. 


3.2.  Infinite  expansion.  The  evolution  of  moments  in  this  case  is  given  by 
the  system  (2.11).  If  velocity  is  constant  then  vj  is  constant  for  each  j,  so  long 
as  the  perturbation  series  can  be  differentiated  termwise.  Thus,  Cy{x,y)  is  also 
constant,  and  Cs^{x,y,t)  is  independent  of  its  second  argument.  Consequently,  is 
identical  to  the  second-order  approximation  to  velocity  variance  and  Cg^{x,  y,  t) 
is  identical  to  Usy{x,t).  Finally,  the  relationship  Ugy  =  UyUg  holds  as  before.  We 
obtain  these  equations  analogous  to  (3.2): 


(3.3) 


dt 


(  So 
e(si) 
(S2) 


\  j 


(  ^^0  /(so)  \ 

^^0  /'(so)e(si) 

{^2)  f{so)  +  cTv  /'(so)o-s  +  /'(so)  (^2) 

+  f'isoWs 

V  (^vf{so)+Vof{so)(Ts  J 


=  0. 


We  have  shown  that  this  system  is  (not  strictly)  hyperbolic  [13].  The  Jacobian 
matrix  in  general  does  not  have  a  full  set  of  linearly  independent  eigenvectors.  This 
degeneracy  leads  to  secular  terms  in  the  solution.  Use  of  the  second-order  mean  as 
the  argument  of  in  §2.3  yields  a  modified  infinite  expansion  that  does  not  have 
this  drawback,  and  is  a  fourth-order  correction  to  the  two-term  expansion  [12]. 


3.3.  Uniqueness,  and  an  additional  analytical  result.  Because  (2.9)  and 
(2.11)  are  nearly  hyperbolic  systems,  one  might  expect  to  extend  uniqueness  meth¬ 
ods  from  the  theory  of  such  systems.  The  viscosity  method  is  an  appealing  way  to 
prove  uniqueness  [25].  Variations  on  this  approach  generally  require  systems  that 
are  genuinely  nonlinear  01  linearly  degenerate  ([25];  see  [4]  for  a  more  recent  result 
and  additional  references).  Neither  the  deterministic  version  of  equation  (2.1)  nor 
the  reduced  systems  above  possess  either  of  these  properties.  Therefore  we  cannot 
apply  existing  uniqueness  arguments  to  (2.9).  A  review  of  the  literature  does  not 
reveal  a  uniqueness  result  general  enough  to  guarantee  uniqueness  for  (2.9).  Thus 
uniqueness  remains  an  open  question;  however,  physical  and  mathematical  argu¬ 
ments  suggest  that  such  results  can  eventually  be  obtained.  For  now,  we  must  be 
satisfied  with 

Conjecture  3.3  (uniqueness).  Moment  equations  (2.9)  or  (2.11)  have  at  most 
one  vanishing-viseosity  solution  for  uniformly  bounded,  measurable  initial  data. 

Now  we  may  obtain  a  more  general  form  of  the  covariance  relationship  stated 
earlier,  directly  from  the  moment  equations. 

Lemma  3.4.  If  there  exists  a  unique  solution  to  (2.9)  with  bounded,  measurable 
initial  data,  then 

(3.4)  Cg  c„  —  Cg„  Cg„ . 

The  proof  in  [13]  uses  the  fact  that  the  same  PDF  is  satisfied  by  both  Cg^  Cg^ 
and  CgC.g.  Again  ag.^  =  aga^  follows,  by  letting  y  ^  x  so  that  Cg  — >■  CTg  and 
Cg.g  Ggy.  It  is  important  to  recognize  that  this  covariance  relationship  follows 
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Figure  1.  Mean  saturation,  linear  flux  with  Gaussian  initial  pro¬ 
file  for  uniform  mean  flow.  Both  the  two-term  (solid)  and  infinite- 
term  (dotted)  moments  exhibit  bimodal  behavior.  The  latter  vio¬ 
lates  physical  bounds  on  saturation. 


from  observation.  That  it  is  also  a  consequence  of  the  moment  equations  (subject 
to  uniqueness)  shows  that  the  MDEs  are  consistent  with  this  intuitive  result. 

4.  Solution 

Solutions  of  nonlinear  hyperbolic  conservation  laws  consist  of  shock  and  rar¬ 
efaction  curves  in  phase  space  [23],  [25].  We  present  results  for  the  linear  advection 
case  first,  then  a  solution  for  nonlinear  advection  with  the  two-term  expansion. 

4.1.  Linear  advection.  Figure  1  shows  mean  saturation  obtained  from  (3.2) 
and  (3.3)  for  linear  fractional  flow  (/(s)  =  s)  with  Gaussian  initial  profile.  This 
linear-flux  case  represents  the  pure-advection  form  of  conservative  solute  transport 
in  single-phase  flow,  with  s  representing  concentration  of  solute.  Both  solutions  are 
bimodal,  and  the  solution  to  (3.3)  violates  physical  bounds  on  saturation. 

This  bimodal  profile  is  a  non-physical  result;  we  do  not  expect  an  initial  local¬ 
ized  pulse  of  solute  to  have  bimodal  behavior  in  the  mean.  Adding  a  linear  diffusive 
effect  will  not  eliminate  bimodality,  as  is  evident  in  the  nonlinear  case  (§4.2),  where 
numerical  and  artificial  diffusion  are  present  in  our  numerical  scheme. 

We  have  used  these  standard  results  for  uniform  mean  flow  for  x  G  [0,  L],  with 
stationary  log  hydraulic  conductivity  ((F)  and  ay  do  not  depend  on  x)  [5],  [10]: 


where  Kq  =  exp((y))  and  J  is  the  negative  mean  head  gradient.  We  used  param¬ 
eters  Kg  J  =  0.5,  (f)  =  0.2,  and  ay  =  0.5. 

Remark  4.1.  Bimodal  mean  concentration  (or  saturation,  in  our  case)  is  noted 
in  [15],  [18],  and  [19].  All  use  methods  to  derive  mean  transport  equations  that  do 
not  involve  second-order  corrections.  Further  details  and  comparison  of  our  work 
to  these  results  can  be  found  in  [12]  and  [13]. 

4.2.  Nonlinear  advection.  Equation  (3.2)  is  in  the  form  (3.1)  with  G  =  0. 
Eor  nonlinear  F(u),  a  solution  consists  of  a  sequence  of  shock  waves,  constant  states, 
and  rarefaction  waves.  We  will  omit  considerable  detail  that  requires  an  appeal  to 
the  elegant  theory  of  hyperbolic  conservation  laws.  Eor  much  more  detail,  see  [13]. 
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(slow  rarefactions) 


Figure  2.  Construction  of  the  slow  rarefaction  Ri  from  (1,0). 


We  refer  to  the  space  with  coordinates  u  =  (mi,  U2)  as  phase  space.  Rarefaction 
waves  are  integral  curves  of  the  vector  fields  defined  by  eigenvectors  of  DF.  Shocks 
are  allowed  if  we  consider  the  weak  form  of  (3.1).  Then  the  Rankine-Hugoniot 
condition  gives  an  expression  for  the  speed  of  the  shock,  and  entropy  conditions 
are  imposed  to  capture  the  vanishing-viscosity  solution.  In  particular,  we  impose 
entropy  conditions  of  Lax  [25]  and  an  extension  due  to  Liu  [21],  [22]. 

We  consider  sure  (deterministic)  initial  saturation  given  by  the  Heaviside  func¬ 
tion  s(x,  0)  =  H[—x\,  so  that  Ug  is  initially  zero.  First  we  simplify  the  notation  of 
(3.2)  as  follows.  Set  ui  =  (s)  and  U2  =  Og,  and  introduce  the  scaling  r  =  {v)  t,  so 
that  dt  =  {v)  dr-  Let  e  =  a^/  {v).  This  ratio  is  consistent  with  the  e  previously 
used  as  the  expansion  parameter.  The  scaled  equations  are 

U2  J  V  ef  +  fU2 


(4.1) 


dr 


f'ul 


=  0. 


Note  that  the  argument  of  the  fcth  derivative  of  /,  is  always  ui.  Let  u  = 
(mi,  M2).  The  Jacobian  matrix  is 


DF{u) 


f  +  ef"u2  +  e/'  +  f''u2  \ 

e/'  +  ru2  f  J  • 


From  §3.1,  its  eigenvalues  are  real  and  may  be  ordered  Ai  <  A2  except  at  u  G 
{(0, 0),  (1, 0)},  where  both  are  zero.  The  associated  eigenvectors  are  r)S;(u),  fc  =  1, 2. 

We  look  for  solutions  to  (4.1)  in  the  phase  space  for  (mi,  M2).  The  points 
u  =  (1,0)  and  u  =  (0,0)  are  endpoints  of  the  solution  to  (4.1).  They  represent  a 
boundary  condition  and  the  initial  condition  within  the  spatial  domain,  respectively. 
We  find  that  a  rarefaction  connects  to  (1,0).  This  curve  Ri  is  shown  in  figure  2 
to  be  an  integral  curve  of  the  eigenvector  associated  with  the  smaller  eigenvalue 
(a  “slow”  rarefaction).  A  shock  must  connect  (0,0).  Using  the  Rankine-Hugoniot 
and  entropy  conditions,  we  construct  the  shock  curve  S2  (see  figure  3). 

The  simplest  connection  between  Ri  and  S2  is  a  single  slow  shock,  denoted 
Si  (see  figure  3).  The  complete  solution  constructed  in  this  manner  is  shown  to 
satisfy  entropy  conditions  in  [13].  The  solution  in  physical  space  is  shown  in  figure 
4,  and  is  compared  to  the  solution  obtained  from  our  numerical  PDF  scheme,  for 
L  =  2,  m  =  1/2,  (v)  =  5/2,  =  5/4  at  a  fixed  time  t  =  0.2. 

Uniqueness  of  the  solution  remains  in  question.  Most  results,  again,  require 
genuine  nonlinearity  or  linear  degeneracy,  and  our  system  does  not  satisfy  these 
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Figure  3.  Construction  of  shocks  connecting  to  rarefaction  Ri 
and  (0, 0).  First,  a  fast  shock  S2  connecting  to  (0, 0)  is  constructed. 
Then  a  slow  shock  is  introduced  to  connect  i?i  and  82- 


Figure  4.  Mean  and  standard  deviation  of  saturation.  The  solid 
curve  is  obtained  from  semi-analytical  construction  of  shocks  and 
rarefaction  curves  in  phase  space;  dashed  curve  is  obtained  from 
numerical  PDF  scheme  applied  directly  to  (4.1). 

conditions.  Liu  extended  existence  and  uniqueness  results  to  more  general  systems 
[21],  [22],  but  the  restrictions  he  places  on  flux  functions  are  not  met  by  Fi  and  F2. 
However,  it  is  clear  in  figure  4  that  this  solution  matches  the  solution  obtained  from 
our  numerical  PDF  scheme,  applied  directly  to  (4.1).  Our  upwind  PDF  scheme  is 
conservative,  and  includes  numerical  and  artificial  diffusion.  The  solution  obtained 
using  this  scheme  is  therefore  an  approximation  to  the  viscous  solution.  The  dif¬ 
fusion  coefficient  is  roughly  three  orders  of  magnitude  smaller  than  the  jumps  in 
solution  values.  Thus  the  numerical  solution  is  near  the  vanishing- viscosity  limit. 
Furthermore,  this  shows  that  linear  diffusion  terms  do  not  eliminate  bimodality  in 
the  solution.  Such  terms  only  smooth  out  mean  saturation  fronts. 

The  saturation  variance  is  supported  primarily  on  an  uncertainty  interval  be¬ 
tween  fronts.  Physically,  the  solution  represents  two  zones  containing  mixtures  of 
the  two  fluid  phases  (for  example,  water  and  oil),  and  a  third  containing  only  the 
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Figure  5.  Buckley-Leverett  saturation  profile  at  t  >  0. 

oil  phase.  In  the  first  zone,  we  have  a  smoothly  varying  mixture  from  the  injection 
boundary  (x  =  0),  where  the  mean  oil  content  tends  to  zero,  down  to  a  constant 
mixture  just  left  of  a  shock.  In  the  second  zone,  we  have  a  constant  mix  of  oil  and 
water.  The  solution  does  not  represent  physical  reality.  Rather  than  two  shock 
waves,  the  true  mean  saturation  is  more  likely  to  have  a  smooth,  diffuse  profile. 
Taken  with  the  profile  of  as,  however,  these  second-order  solutions  provide  some 
insight  into  the  propagation  of  uncertainty  in  two-phase  flow. 

In  the  limit  a^  0,  the  mean  saturation  tends  to  the  classic  Buckley-Leverett 
profile  shown  in  figure  5.  In  fact,  the  construction  of  the  rarefaction  and  shock 
profile  is  analogous  to  the  solution  of  the  scalar  deterministic  saturation  equation. 
For  that  solution,  an  initial  rarefaction  is  followed  forward  in  x,  up  to  a  point  where 
the  shock  speed  matches  the  characteristic  speed. 

5.  Conclusions 

For  spatial  resolution  of  uncertainty,  these  bimodal  results  suggest  that  second- 
order  Eulerian  MDEs  may  be  inappropriate  for  1-D  immiscible  flow.  Unlike  1-D,  v 
has  finite  correlation  length  in  2-D,  and  macrodispersion  models  are  based  on  cor¬ 
relations  in  5v.  In  results  for  analogous  MDEs  for  passive  2-D  solute  transport  with 
diffusion,  bimodality  was  not  observed  [11];  one  might  expect  macrodispersion  to 
lead  to  a  similar  result  for  immiscible  flow.  Nevertheless,  in  a  forthcoming  submis¬ 
sion  we  show  that  somewhat  mitigated  bimodality  does  persist  in  2-D,  even  with 
diffusion  terms  [14].  More  positively,  for  spatial  averages  such  as  oil-production 
curves,  good  matches  to  Monte  Carlo  simulations  are  found. 
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